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Abstract. We use a generic model for type-I excitability (known as the SNIPER or SNIC model) to de- 
scribe the local dynamics of nodes within a network in the presence of non-zero coupling delays. Utilising 
the method of the Master Stability Function, we investigate the stability of the zero-lag synchronised 
dynamics of the network nodes and its dependence on the two coupling parameters, namely the coupling 
strength and delay time. Unlike in the FitzHugh-Nagumo model (a model for type-II excitability), there 
are parameter ranges where the stability of synchronisation depends on the coupling strength and delay 
time. One important implication of these results is that there exist complex networks for which the adding 
of inhibitory links in a small-world fashion may not only lead to a loss of stable synchronisation, but may 
also restabilise synchronisation or introduce multiple transitions between synchronisation and desynchro- 
nisation. To underline the scope of our results, we show using the Stuart-Landau model that such multiple 
transitions do not only occur in excitable systems, but also in oscillatory ones. 



1 Introduction 

Studies on complex networks have gained much atten- 
tion in recent times from various fields of research [11213] . 
In particular, the dynamics on networks of nonlinear ele- 
ments (or nodes) has become a very active area 
Here, we consider the dynamics of neurons as a network of 
nonlinear excitable nodes. While the model we consider is 
generic for excitability, the local dynamics being modelled 
in this context is the "firing" mechanism of neurons (i.e. 
the rise and fall of the electrical membrane potential). This 
mechanism is the basis for neural coding and information 
transfer or cell-to-cell communication [15 . The type of 
network we primarily focus on in our investigations is the 
small- world (SW) network [16] . which has a short average 
path distance between nodes, as well as a large degree of 
clustering (in other words, many triangles in the network 
structure). These properties are found in many kinds of 
real- world structures, such as the collaboration between 
film actors, power grids, the World Wide Web and social 
relationships |16|17|1|2] , In particular, large-scale cortical 
networks also show these properties |18|19j . The brain has 
an architecture enabling both efficient global and local 
communication between neurons [20] . which is captured 
well by the SW model. 



The interest in the synchronisation of neurons emerges 
from the role it plays in information processing and per- 
ception, as well as in epilepsy and similar diseases |21|22|23|24|25j . 
While there are several types of synchronisation, the fo- 
cus here will be on zero-lag synchronisation. This means 
tfr flit ^U p tly T J ^ rn ^ rs °f a U nodes in the network are not only 
fit also temporally in-phase. 
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Inhibition plays an important role in the nervous sys- 
tem [26] . Here, when constructing a network we begin 
with a regular ring network of excitatory links and, as 
in Ref. [8], we add long-range inhibitory links into the 
network structure. This creates a SW network of the form 
proposed in Ref. [27] . The introduction of inhibitory links 
and an increasing of their density in the network was 
shown in Ref. [8 to result in a loss of stable synchro- 
nisation. 

In this paper, we consider delay-coupled networks of 
type-I excitability, in contrast to Ref. [8] where type-II 
was considered using FitzHugh-Nagumo local dynamics, 
and find qualitative differences for a particular range of 
coupling parameters. The most interesting implication of 
this is that for some network topologies with this type 
of excitability increasing the density of inhibitory links 
may not just lead to one transition from stable to unsta- 
ble synchronisation, but may lead to multiple transitions. 
This means that increasing the density of inhibitory links 
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in the network may stabilise synchronisation, rather than 
destabilise it. 

In the following section, we introduce the model and 
the network topologies considered in the present study. In 
Sec. [3] we use the master stability function to investigate 
the stability of arbitrary synchronised networks with given 
coupling parameters (i.e. the coupling strength and the 
length of delay between coupled nodes). The master sta- 
bility function is calculated in Sec. [4] for networks coupled 
within a range of small delay time and coupling strength, 
which delivers the main results of the paper - namely, the 
existence of synchronised states that have different sta- 
bility conditions compared to coupling with larger delay 
times. The implications these results may have for specific 
complex networks are discussed in Sec. [5] In Sec. |6j we 
study networks of coupled Stuart-Landau oscillators and 
demonstrate that multiple transitions between synchroni- 
sation and desynchronisation can also appear in networks 
of oscillatory nodes and are not limited to excitable sys- 
tems. 



2 Model 

In order to model excitability, the system must have a 
rest state, which corresponds to a stable fixed point. Small 
perturbations from the rest state allow for the creation of 
a large excursion in the phase space, which corresponds 
to the system entering an excited state before returning 
to the rest state. In the context of neurodynamics, this is 
the firing state of the neuron [28 . 

Neurons can exhibit different excitability properties. 
In 1948, Hodgkin classified two types of neural excitability 
[29] . Both types can be ordered depending on the type of 
bifurcations that occurs when the bifurcation parameter 
is changed such that the system makes a transition from 
the excitable to the oscillatory regime: 

Type-I neurons can generate action potentials of arbi- 
trarily low frequency. This kind of behavior occurs near 
a saddle-node infinite period (SNIPER) bifurcation, also 
known as the SNIC bifurcation (saddle-node bifurcation 
on invariant cycle). The arbitrarily low frequency coin- 
cides with the period of the limit cycle going to infinity 
as the bifurcation parameter approaches a critical value, 
where the bifurcation occurs. 

Type-II neurons are associated with a supercritical 
Hopf bifurcation. The frequencies of the action potentials 
lie within a certain non-zero range, while the amplitude 
of the limit cycle approaches zero with the bifurcation pa- 
rameter. 

As our model for type-I excitability we consider a generic 
normal- form of a SNIPER bifurcation [30 31 . This model 
is mathematically represented by two first-order differen- 
tial equations: 

fM- f x \ - fx(l~x 2 -y 2 )^y(x-b)\ m 
W ~ W 1 ~ x2 - V 2 ) -x(x-b)J' [ j 



In polar coordinates (x = r cos (p,y = r sin (p) Eq. ([!]) reads 

r = r(l-r 2 ) (2) 
(f = b — r cos (p. (3) 

b > is the bifurcation parameter. This bifurcation pa- 
rameter influences the type of dynamics and determines 
where in the (x, ?/)-plane the fixed points are located, as 
discussed below. 

An important feature that must be taken into account 
in the mathematical model is that the effect one neuron 
has on another is not instantaneous. The signals sent be- 
tween neurons have a finite propagation speed and take 
a certain time to arrive at the destination neuron, giving 
rise to a delay time r. 

The dynamics of a network of N elements (labelled 
i = 1, N) is written: 

N 

x, = f (xi) + a G %J K(xj(t - r) - x,(t)), (4) 

where f(x^) is the local dynamics as described by Eq. 
for each element x^ = (x^yi). G^ determines the matrix 
G for the network structure, showing which elements are 
coupling together, and H is the coupling function. H is 
taken to be the 2x2 identity matrix; this means that 
the x- variable at time t is coupled with the x- variable at 
time t — r, the y- variable at time t is coupled with the 
2/- variable at time t — r, but there is no cross-coupling 
between the x- and variable. The coupling parameters, 
which are identical for all connections, are the coupling 
strength a and delay time r. 

For synchronous dynamics, the 2 (TV — 1) constraints 
x, s = Xi = x 2 = • • • = xat define a 2-dimensional synchro- 
nisation manifold (SM) within the 2iV-dimensional phase 
space and the coupled system is reduced to an effective 
single system with feedback: 

x s = f (x s ) + (iH(x s (t - r) - x a (t)), (5) 

with unity row sum . Gij = 1, so that the nodes all 
receive the same level of input while they are synchronised. 
Any non-unity but constant row sum can be rescaled using 
the coupling strength a. 

With the bifurcation parameter b < 1 there exists an 
unstable focus at the origin, as well as a saddle point 
and a stable node situated on the unit circle at {x^yi) = 
(6, y/l — b 2 ) and (6, — \A — fr 2 ), respectively. At b = 1 the 
saddle point and stable node collide, so that for b > 1 
a limit cycle exists along the unit circle. In the case of 
b < 1 (excitable regime), a perturbation which pushes the 
system from the stable node beyond the saddle point can 
result in a single oscillations along the unit circle. Delayed 
coupling can induce a homoclinic bifurcation, such that a 
limit cycle is produced that bypasses the saddle point and 
stable node [31 32 . Here, we will consider the excitable 
regime with b = 0.95 and investigate for which topologies 
the synchronised dynamics of complex networks is stable. 

In this paper, two kinds of complex networks are con- 
sidered: SW networks and Erdos-Renyi random networks 
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[33] . We construct the SW networks as a variation to the 
method proposed in Ref. [IB] , introduced in Refs. |34|27j : 
(i) Each of the TV nodes in a one-dimensional ring is given 
excitatory links to its k nearest neighbours on each side. 
Note that in terms of the matrix G, an excitatory link be- 
tween the i th and j th node means that Gij > 0, while for 
an inhibitory link Gij < 0. (ii) For each of the kN links 
we add with a probability of p another inhibitory link 
with strength -1 between two randomly chosen nodes (i.e. 
on average pkN randomly distributed inhibitory links), 
(iii) Self-coupling and multiple links between the same two 
nodes are not allowed, (iv) Finally, the entries in each row 
of G are normalised to ensure a unity row sum. If a row 
sum is equal to zero, then the network realisation is dis- 
carded. 

The random networks are constructed with kN exci- 
tatory links between randomly chosen nodes, as well as 
pkN randomly distributed inhibitory links, so that these 
may be compared to a SW network with the same param- 
eters p, k and N. Because the construction of both kinds 
of complex networks involves random processes, there will 
be many possible realisations of either kind of network 
with a given set of parameters. 

In the following, we determine stability of the synchro- 
nised dynamics and compare the stability of these different 
kinds of networks. 
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Fig. 1. (Colour online) Master stability function A for 
coupling parameters cr = 0.3 and (a) r = 10, (b) r = 7, 
(c) r = 6.5, and (d) r = 6. a and j3 are the real and imagi- 
nary parts of the scaled eigenvalues of the coupling matrix, 
cr^/e, respectively, and A is the largest Lyapunov exponent. 
White dots in panel (a) represent the scaled eigenvalues of 
a unidirectional ring of 11 nodes (crz/ l5 . . . , av\\). b = 0.95. 



3 Stability of Synchronisation 

The Master Stability Function (MSF) is a tool that repre- 
sents a measure of the stability of a synchronised state for 
given delay-coupling parameters in terms of the topology 
of an arbitrary network [35]. The master stability equa- 
tion, given by 

6i(t) = [Df (x s ) - o-H](Jx(t) + (a + z/?)H£x(* - r), (6) 

is found by linearising Eq. Q around Eq. ([5| and can be 
used to calculate the largest Lyapunov exponent A(a, /?, cr, r) 
called the MSF. is the perturbation of x away from the 
SM (Le. x = x s + (5x) and Df(x s ) is the Jacobian matrix of 
Eq. Q evaluated on the SM. Important for this approach 
is that, whereas one can calculate Lyapunov exponents for 
a specific network topology using the eigenvalues of the 
matrix G, one considers here a continuous complex pa- 
rameter a + if3 representing the complex plane of possible 
eigenvalues scaled by the coupling strength cr (i.e. a + i/3 
is a continuous parametrisation of ovj, where Vj are the 
eigenvalues of G, j = 1, N). Thus, one can calculate the 
Lyapunov exponents for a region of the (a, /3)-plane which 
gives sub-regions of stability, where A < 0, and instability, 
where A > 0. It is then easy to compare the synchronous 
stability of various networks by simply observing whether 
any of their eigenvalues fall inside an unstable region of 
the (a, /3)-plane. If all the eigenvalues lie within stable re- 
gions, then perturbations away from the SM will decay 
exponentially. 

Because of the unity row sum condition, G always has 
an eigenvalue v\ — 1. This longitudinal eigenvalue (all 



others are called transversal) corresponds to perturbations 
within the SM and A{av\^ 0, cr, r) is always zero because 
we are looking at periodic dynamics. As such, it is only the 
transversal eigenvalues that are important for determining 
the stability of synchronisation. 

Figure^a) shows the MSF of the SNIPER model with 
coupling parameters a = 0.3 and r = 10. The white dots 
represent the eigenvalues of the matrix G for a unidirec- 
tional ring of 11 nodes. All eigenvalues are within the sta- 
ble region of the MSF, thus the synchronisation of all 11 
nodes coupled with these parameters is stable. 
} According to Ref. [7], for r in the order of the system's 
characteristic time scale (in this case, the period of the 
oscillations) and above, the MSF will always tend towards 
a rotational symmetry. Examples such as the one above in 
Fig. HI a) confirm these general findings for the SNIPER 
model. When calculating MSFs for a large fixed r while 
varying cr it becomes evident that the size of the stable 
region is scaled by cr, so that the stable region can be 
estimated very well by the circle £((0,0), cr) (that is, a 
circle centred at the origin with a radius cr). Large values 
of r, however, despite influencing the Lyapunov exponents 
quantitatively, do not affect the shape of the stable region. 

This rotational symmetry of the MSF was also found 
for type-II neurons, modelled as FitzHugh-Nagumo oscil- 
lators [8]. In this case a and r do not affect whether the 
eigenvalues fit into the stable region of the MSF, so that 
only the topology of a network is important for the sta- 
bility of its synchronisation. Furthermore, because Gersh- 
gorin's circle theorem [36 guarantees that the eigenvalues 
of a network's coupling matrix with no self-coupling and 
purely excitatory coupling (i.e. Gu = and Gij > 0, 1 < 
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Fig. 2. (Colour online) Fraction of desynchronised net- 
works f(p) vs. probability of inhibitory links p for 100 
realisations of SW networks of N = 100 and k = 24 for 
a = 0.3, r = 10, 6 = 0.95. 

N) lie within the unit circle on the complex plane, 
the synchronisation of such a network will always be sta- 
ble. Finally, if additional inhibitory links are introduced 
with probability p to construct a small-world network as 
described in Sec. [2j phase transitions from stable to un- 
stable synchronisation are found with increasing probabil- 
ity of inhibition p [8 . All these results also apply for the 
SNIPER neurons with sufficiently large r. 

As an example, Fig. [2] shows the transition from stable 
synchronisation, where the fraction / of network realisa- 
tions with unstable synchronisation is zero, to unstable 
synchronisation, where / = 1, as the probability of in- 
hibition p is increased. This example is for the coupling 
parameters a = 0.3 and r = 10, but this transition to 
inhibition-induced desynchronization will look the same 
for all sufficiently large values of r. A thorough investiga- 
tion of the effects of additional inhibitory links when r is 
small will be presented in Sec. [5| 



4 Small delay times 

If the delay time is not large enough, the rotational sym- 
metry of the MSF no longer holds. In fact, while reducing 
r one can witness how the rotational symmetry begins to 
break down. This is depicted in Fig. [I] For a constant cou- 
pling strength of a = 0.3, the MSFs are numerically calcu- 
lated for decreasing delay times. While at r = 10 the MSF 
still has its circular form (Fig. [TJa)), when decreasing r, 
the stable (i.e. dark blue/green) region of the MSF begins 
to show signs of deformation. By r = 7 (see Fig.JIJb)) the 
stable region is clearly larger than the unit circle scaled 
by a — 0.3 and has definitely lost its rotational symmetry. 
By r = 6.5 (see Fig. [TJc)) the stable region has already 
split into two disconnected regions. Letting r decrease fur- 
ther, the stable regions become increasingly smaller (see 
Fig.JIJd)). Note that the delay-induced limit cycle coexists 
alongside the stable fixed point and whether it is reached 
or not is therefore dependent on initial conditions. For r 
less than about 4 (not shown here), the coupling is no 



longer able to induce the homoclinic bifurcation that cre- 
ates the limit cycle at all (as discussed in Sec. [2]), in other 
words, the neurons no longer oscillate. Instead, all solu- 
tions approach the stable fixed point. 

It is now obvious that, in this regime of small delay, 
small changes in r can have a great impact on the stabil- 
ity. Despite the seemingly sudden change in the MSFs in 
Fig. [I] between r = 7 and r = 6.5, the evolution of the 
boundary of stability is in fact a continuous one (except 
at a critical value r c , which will be discussed below). This 
becomes clear after plotting the MSF versus the real part 
a of the eigenvalue (with f3 = 0) for varying r. Taking this 
one slice of the eigenvalue plane gives a good indication 
of the growth and decay of the stable region in the MSF 
while changing r. Figure [3] shows the MSF as a function 
of the real part a (/3 = 0) and the delay time r, with 
fixed coupling strengths of (a) 0.25, (b) 0.30 and (c) 0.35. 
This produces an MSF on an (a, r) -plane, which can be 
used for network topologies with a coupling matrix that 
yields real eigenvalues, e.g. symmetric matrices for undi- 
rected networks. In the following we restrict ourselves to 
undirected networks. 

Based on the exemplary coupling strengths in Fig. [3] 
it is clear that the stability depends on both a and r. In 
all examples the vertical boundary line at a = <r, cor- 
responding to the longitudinal eigenvalue (i.e. vq = 1, 
where A = 0), can be easily identified. It separates regimes 
of stable and unstable synchronisation. Another common 
characteristic is that the r-dependent MSFs have a crit- 
ical delay time r c at which the stable region splits into 
two separate, disconnected regions. For values above r c 
the stable region is found to the left of the longitudinal 
eigenvalue; whereas for values below r c there may be one 
stable region to the right of the longitudinal eigenvalue 
and one to the left. r c seems to be an important value, 
because it marks the most significant r-dependent qual- 
itative change in the MSF. Ultimately, the MSF can be 
divided into three regimes of r: (i) two or more smaller 
separate regions of stability exist when r < r c ; (ii) one 
large region of stability exists when r > r c ; (iii) the MSF 
has one rotationally symmetric region of stability in the 
limit of r — » oo (which holds already in good approxima- 
tion if r is of the order of the intrinsic oscillation period 
or larger). 

Clearly, the topological criteria for stable synchroni- 
sation will be very different depending on whether the 
network is coupled with r < r c or r > r c . Observing the 
change in the limit cycles while varying r offers the ex- 
planation for the qualitative change either side of r c in 
the MSF. For smaller r values, the trajectory of the limit 
cycle is sensitive to changes in the coupling parameters, 
whereas, for larger r values, the trajectory changes very 
gradually and approaches an elliptical form. The most pre- 
dominant feature of the changing limit cycle is a large 
"kink" for values below r c . Figure pla) shows such a kink 
in the limit cycle induced by the coupling parameters 
<j = 0.35 and r = 5.35: upon bypassing the stable node 
the limit cycle deviates a little towards the unstable focus, 
located at the origin, before suddenly changing its trajec- 
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Fig. 3. (Colour online) MSF A for a fixed coupling 
strength (a) a = 0.25, (b) a = 0.3 and (c) a = 0.35 in 
the plane of the real part a (/? = 0) and the delay time r. 
The horizontal red lines show the position of the critical 
delay time r c . b = 0.95. 



tory back towards the saddle point. The kink has a large 
effect on the dynamics of synchronised nodes and is an im- 
portant source of instability. Simulations have shown that 
it is at the kink that synchronised nodes begin to diverge 
from one another following any small perturbations. For 
larger r the kink diminishes and disappears (see Fig.Qb)). 
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Fig. 4. (Colour online) The position in phase space of 
one element representing the synchronised dynamics of a 
network (large green diamond) is shown relative to its 
delayed position (small green diamond). The black circle 
and the white circle represent the stable and saddle point, 
respectively. For small delay times such as (a) r = 5.35 
a kink may be the result of the relative positioning. For 
larger delay times such as (b) r = 6 this does not occur. 
The movement of the phase point has a local component 
(dark green arrow) and a coupling component (light green 
arrow). The coupling component is directed towards the 
phase point at time t — r (small green circle), a = 0.35, 
6 = 0.95. 



Let us explain the cause for the kink in greater de- 
tail. The large green diamond in Fig. [4ja) is the position 
of a node that represents the dynamics of a synchronised 
network in the phase space. The small green diamond is 
its delayed phase point, to which it is coupled, and rep- 
resents the delayed dynamics of the same synchronised 
network (i.e. the small green diamond at time t is the 
large green diamond at time t — r). A node coupled with 
a delayed node is attracted towards the direction of this 
delayed node (cf. Fig. [4^b)). As the node passes between 
the stable fixed point and the saddle point, it is slowed 
down, allowing its delayed node to "catch up" to it. For 
smaller r this has the effect that the node is attracted in 
a clockwise direction towards its delayed node and gets 
slowed down (cf. Fig.Qa)). This is the cause for the kink. 
Whereas for larger r, the (small green) delayed node is 
always more than half a rotation behind, meaning that 
the (large green) node is attracted towards it in an anti- 
clockwise direction. This gives the limit cycle a form that 
is roughly circular as seen in Fig. [4^b). 

As long as the kink is present in the limit cycle, a syn- 
chronised network will generally have very different topo- 
logical conditions in order to be stable. Even when the 
limit cycle has a very predominant (sharp) kink, it will 
often be possible to find a class of networks that will still 
show stable synchronisation. 

Note that in Fig. [3] one can see how the most stable 
part of the stable region (i.e. where A is at its smallest) 
shifts away from a = as r becomes smaller. Simulations 
have shown that this happens while the angle of the kink 
increases (making the kink "sharper"). A possible inter- 
pretation could be as follows. Each potential eigenvalue 
of the coupling matrix in a MSF is associated with one 
or more eigenvectors in the A^-dimensional space. These 
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eigenvectors are all orthogonal to the vector (1, 1, 1, 1) 
that represents motion in the SM. This would suggest that 
the increasing angle of the kink corresponds to a shifting 
of the most stable direction of perturbation transverse to 
the SM. 



5 Implications for complex networks with a 
small delay coupling time 

An obvious observation is that networks with purely ex- 
citatory coupling, which are always stable for large delay 
times as mentioned at the end of Sec.[3j may not show sta- 
ble synchronisation for r < r c . It was explained in Sec. [3] 
that increasing the probability p of inhibitory links in the 
network can be a factor leading to unstable synchronisa- 
tion. This occurs because a larger probability of inhibition 
can push a part of the eigenvalue spectrum of any network 
beyond the longitudinal eigenvalue at a = a (because Ger- 
shgorin's circle theorem no longer holds) and into the un- 
stable region of the MSF. Now, for smaller r, there may be 
a pocket of stability to the right of the longitudinal eigen- 
value, so that increasing inhibition can make the otherwise 
unstable synchronisation of a network stable. 

This means that the transitions between stable and 
unstable synchronisation as a function of the probability 
of inhibition p discussed in Ref. [8] are no longer valid 
when r is small. The transitions are now sensitive to the 
coupling parameters, not just the network parameters N 
and k. Due to the multiple regions of stability, eigenvalues 
may wander in and out of stable regions, while increasing 
the probability of inhibition p. For large r, increasing p 
in the SW network model only results in one transition 
where the fraction of desynchronising networks (i.e. net- 
works with unstable synchronisation) f(p) switches from 
to 1 (cf. Fig.[2|. For small r, it is possible that f(p) jumps 
back to 0, before increasing again to 1. This will occur 
if there is a separate region of stability to the right side 
in the MSF that is large enough that all the eigenvalues 
lying over there can fit inside. 

The observation of multiple transitions between stable 
and unstable synchronisation can be explained by looking 
at the eigenvalue spectra for SW networks for various p 
values. As discussed above in relation to the MSF method, 
each network topology has a coupling matrix G, the eigen- 
values from which can be used to determine the stability 
of the network's synchronised state. Because the "short- 
cuts" in the SW network are randomly introduced, each 
network with specific TV, k and p values may have many re- 
alisations. By calculating the eigenvalues for a large num- 
ber of realisations of the SW network with certain param- 
eters (i.e. given N, k and p), one can find the bounds for 
the eigenvalue spectra. 

Figure [5] displays the superimposed eigenvalue spectra 
of 500 realisations for SW networks with parameters N — 
200 and k = 20 (regular ring with excitatory coupling of 
k nearest neighbours on either side) in dependence on the 
probability of inhibition p. The longitudinal eigenvalues 
located at ui = 1 have been removed. One can see that 
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v 

Fig. 5. Spectra of eigenvalues v of undirected small world 
networks with N = 200 and k = 20 and varying probabil- 
ity p of inhibitory links. 500 realisations of the eigenvalue 
spectrum are plotted for each value of p. 

the bounds for possible eigenvalues shift depending on the 
value of p. One can also see how the spectra begin to 
increasingly resemble the semicircular distribution pQ of a 
random network for larger p values, where the networks 
have lost their SW properties. 

The eigenvalue spectra bridge the gap between obser- 
vations of the MSF (i.e. the dependence of the stability 
of synchronisation on the eigenvalues) and what is seen in 
these transitions (i.e. the dependence of stability on the 
network topology). Increasing p allows isolated eigenval- 
ues to increase in value and, in terms of the MSF, shift 
their locations further to the right in the (a, /3)-plane. This 
is visualised in Fig. [6] Figures [6^b), (d) and (f) show the 
eigenvalue spectra for SW networks with N = 200 ele- 
ments and k = 40, 20, and 50, respectively. 

Figures^ a), (c) and (e) show the corresponding frac- 
tion of desynchronised networks / in dependence on the 
probability of additional inhibitory links p. Consider, for 
instance, SW networks with parameters N = 200 and 
k = 40. In Fig. |6ja) f(p) is shown for the exemplary 
coupling parameters a = 0.3 and r = 6.5 with a cor- 
responding plot of eigenvalues in Fig. [6lb). Note that the 
grey (green) shaded regions represent the stable regions of 
the line ft = on the (a, /3)-plane of the MSF for a = 0.3 
and r = 6.5 (cf. Fig.[TJc)). For coupling parameters where 
this region of stability is not large enough, / may briefly 
dip below 1 without decreasing to 0, because only some 
realisations may have eigenvalues that fit inside the stable 
region. In Fig. |6jc) where k = 20, / dips down to 0.14, 
because at most 14% of the realisations show stable syn- 
chronisation (i.e. all the eigenvalues are inside the stable 
region). If the distance between the larger isolated eigen- 
values matches the distance between stable regions (which 
is almost the case in Fig. [6jd)), then the transition curve 
only just touches the / = axis at some value of p be- 
fore increasing back to / = 1. Figure |6Fe) is an example 
where a further transition is possible, because not only 
do all eigenvalues begin in stable regions at p = 0, but 
there is another regime of p where all eigenvalues fit in- 
side stable regions. Note that the length of the transition 
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Fig. 6. (Colour online) Fraction of desynchronized networks / in dependence of the probability of additional inhibitory 
links p for 500 realisations of networks of N = 200 with (a) k = 40, (c) k = 20, and (e) k — 50 with <j — 0.3 and 
r = 6.5. Corresponding eigenvalue spectra for (b) k = 40, (d) k = 20, and (f) k = 50 with 100 realisations for each p 
value. Here the green shaded regions represent the stable regions of the real part of the MSF. b — 0.95. 
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Fig. 7. Histograms of eigenvalue spectra for 1000 realisa- 
tions of (a) SW networks and (b) random networks with 
the same number of nodes, the same number of excita- 
tory links and the same probability of inhibitory links. 
TV = 100, k = 10, p = 0.2 and number of bins is 1000. 



from / = to / = 1 is actually a measure of the variance 
of the isolated eigenvalues for an ensemble of realisations. 
For example, the variance of the isolated eigenvalues de- 
creases as N is increased, so that the length of transition 
will be shorter in larger networks, and the transition be- 
comes sharper. 

When constructing a histogram of the eigenspectra for 
a particular value of p - see Fig[7^a) - the larger isolated 
eigenvalues seen to the right, e.g. in Fig. [6^d), result in 
small peaks of eigenvalues. For each realisation of many 
SW networks there are isolated eigenvalues that are larger 
than most other eigenvalues in the spectrum. The small 
peaks of eigenvalues result from perturbations in isolated 
eigenvalues. The spectrum of a regular ring network (i.e. 



a SW network with p = 0) can be found analytically using 
the graph's symmetry operations [37] and is given by 



1 k 

h 



cos I 2irj 



I 



1 / cos (kirjj) sin ((fc + 1)tt^) 



sin (n 



(7) 



N J 



where / = 1,...,7V — 1. Figure [8] shows this eigenspectrum 
for the p = case. As p is increased these eigenvalues 
will be slightly perturbed (in a random manner) by the 
changing network structure, of which many realisations 
are possible. In Fig. [5] one can see the isolated eigenvalues 
to the right-hand side of the plot that eventually evolve 
into the smaller side peaks of eigenvalues in histograms for 
larger p (see Fig [7^ a)). These peaks show the distribution 
of the eigenvalue under the influence of the random nature 
of the SW "short-cuts" creation process - which is why 
the smaller eigenvalue peaks can be approximated by a 
normal distribution^] This is furthermore the reason why 
the desynchronisation transitions observed in Ref. [8] and 
Figs.[6ja), (c) and (e) look like the cumulative (integrated) 
distribution function. They reflect how increasing p brings 



1 Actually, because the multiplicity of the above mentioned 
isolated eigenvalues at p = is 2, the smaller eigenvalue peaks 
are two normal distributions that, at least for small p values, 
overlap each other to a large extent 
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Fig. 8. The eigenvalue spectrum for a regular network 
with N = 100 and k = 10. v n is the value of the n th 
eigenvalue. 



the area of the small eigenvalue peak accumulatively into 
the unstable region of the MSF. 

To explain why the peaks wander with increasing p, 
consider not normalising the rows of matrix G, so that the 
row sums are not necessarily equal to 1. Then the location 
of the longitudinal eigenvalue decreases with p, because 
it is equal to the average row sum of G which is equal 
to 2fc(l — p). The locations of the other peaks increase 
slightly to maintain the eigenvalue sum of zero (a result 
of the trace of G being zero, since there is no self-feedback 
coupling), although for large N this has only a small ef- 
fect. Thus, scaling the eigenvalues so that the longitudinal 
eigenvalue is always at 1 means that the transversal eigen- 
values are multiplied by 2 k(i-p) ? so ^at they appear to 
increase with p, as seen in Figs.[6jb), (d) and (f). 

The case of random Erdos-Renyi networks [33] is dif- 
ferent. At p = there are no larger gaps between the 
eigenvalues, such as those for SW networks at p = 0. If a 
random network is fully connected (i.e. is one single com- 
ponent) then the distribution of its eigenvalues at p = 
are confined to a semi-circle, given by the semicircular 
distribution 



p{y) = <^ 2<KNq{l-q) 

lo 



if h < 2 
otherwise, 



y/Nqjl-q 
q(N-l) 



(8) 



where q is the probability of excitatory links [1 j. While 
there are initially no gaps in the eigenvalue distribution at 
p = 0, adding inhibitory links to the network (i.e. increas- 
ing p) does not lead to the creation of any gaps like those 
seen for SW networks. To illustrate this point, Fig. [7^b) 
shows the eigenvalue distribution for many realisations of 
Erdos-Renyi random networks with the same number of 
nodes, as well as the same number of excitatory links and 
the same probability of inhibitory links. As a result, mul- 
tiple transitions between stable and unstable synchronisa- 
tion will not occur. 

One can, however, use Eq. ([8| to find parameters N 
and q such that the eigenvalue distribution is confined to 
a stable region of the MSF at p = 0, as long as there is 



a stable region centred around v — (for example, see 
Fig. [TJc)). This allows for a desynchronisation transition 
with increasing p, such as in Fig. |2j 



6 Small-world networks of Stuart-Landau 
oscillators 

The appearance of multiple transitions between synchro- 
nisation and desynchronisation is based on the interplay of 
the peculiar spectral properties of the small-world topol- 
ogy and the unconnected stable regions that show up in 
the MSF for the excitable SNIPER model for small delay 
times. 

In this Section we show that such unconnected stable 
regions may also appear in the MSF of a purely oscillatory 
system. We consider the Stuart-Landau oscillator, which is 
the normal form of any oscillatory system near a supercrit- 
ical Hopf bifurcation. The dynamics of the Stuart-Landau 
oscillator is governed by a complex variable z G C, such 
that we have = zi in Eq. Q. The local dynamics of the 
Stuart-Landau oscillator is given by 



f(z) = (\ + iuj-\z\ 2 )z. 



(9) 



For A > 0, the uncoupled system exhibits self- sustained 
limit cycle oscillations with radius ro = \/\ and frequency 
uj. The coupled system shows - depending on coupling 
parameters and topology - a variety of synchronised and 
desynchronised solutions, including amplitude death and 
cluster synchronisation |38|39j . 

Here we focus on complete (zero-lag) synchronization. 
Fig. |9ja) shows the MSF for sets of parameters that allow 
for unconnected stable regions in the complex plane of 
the eigenvalues (a,/?). As shown in Fig. [9jb) a connected 
stability region occurs for integer multiples of r = 2tt/uj = 
27r, while unconnected stability regions - which can lead 
to multiple phase transitions in SW networks - occur in 
between. 

For the parameters used in Fig. [9ja), Fig. [9^c) shows 
the fraction of desynchronised networks for the modified 
small- world model as in Sec. |5j Figure ^d) shows the 
eigenspectra for networks of N = 200 elements, starting 
from a one-dimensional ring with excitatory coupling to 
k = 40 neighbours to either side. Increasing the probabil- 
ity of additional inhibitory links p broadens the spectra 
and leads to a second regime of stable synchronisation 
matching the transitions in Fig. [9jc). 



7 Conclusion 

We have investigated transitions between synchronisation 
and desynchronisation in complex networks of delay-coupled 
excitable elements of type I, induced by varying the bal- 
ance between excitatory and inhibitory couplings in a small- 
world topology. In our analysis we have used the master 
stability function approach. For large delay times it seems 
that both type-I neurons, as modelled here, and type-II 
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Fig. 9. (Colour online) Network of Stuart-Landau os- 
cillators according to Eqs. Q and (a) MSF in the 
(a, /?) plane for r = 37r/2. (bjMSF in the (a,r) plane for 
P = 0. (c) Fraction of desynchronized networks / in de- 
pendence of the probability of additional inhibitory links 
p for 500 realisations of networks of N = 200 and k — 40, 
r = 37r/2. (d) Corresponding eigenvalue spectra, where 
the green shaded areas represent stable regions. Other pa- 
rameters: A = 0.1, u) = 1, a = 0.08. 



neurons, as modelled in Ref. [8], must fulfill similar topo- 
logical conditions in the network to allow for a stable syn- 
chronised state. This is different when considering small 
delay times. For a range of small coupling strengths and 
small delay times we have found novel multiple transitions 
between synchronisation and desynchronisation, when the 
fraction of inhibitory links is increased. Unlike the Erdos- 
Renyi random network, a small world model for complex 
networks with regular excitatory couplings and random 
inhibitory shortcuts has eigenvalue spectra with gaps be- 
tween the larger eigenvalues, so that histograms of many 
realisations reveal isolated peaks of possible eigenvalues. It 
was shown that, because of this, synchronised small world 
networks can go through multiple transitions of stability 
in dependence on the probability of inhibitory short-cuts. 
Our results are valid beyond the special model of type-I 
excitability, as we have demonstrated using the Stuart- 
Landau oscillator. 

This work was supported by the DFG in the framework 
of the SFB 910. PH acknowledges support by the BMBF 
(grant no. 01GQ1001B). 
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